// 2022-7-14 将不同投影的数据拼接到同一个输出数据上，输出数据为WGS84

#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include "wstringutils.h"
#include "wGdalRaster.h"
#define WGDALRASTER_H_IMPLEMENTATION
#include "wGdalRaster.h"
#include <iostream>
#include "gdal_priv.h"
#include "ogr_spatialref.h"
#include <memory>
#include "wGdalRasterLiner.h"


using namespace std;

struct LLValue{
	double x,y,ndvi,vza ;
};

 
bool ReadFileList(string listfilename ,vector<string>& filenameVec)
{
	std::ifstream infs(listfilename.c_str());
	if( infs.good()==false ) return false ;
	std::string line;
	while (std::getline(infs, line))
	{
		if( line.length() > 2 )
		{
			filenameVec.push_back( line ) ;
		}
	}
	return true ;
}

//使用两次填充，第一次in->out 记下四角范围 第二次从四角范围出发进行 out->in 填充
bool FillOneFileIntoOutput( 
	wGdalRasterLiner* outras ,
	string infile1 ,
	double filldata,
	double vmin,double vmax,
	string& error
	)
{
	shared_ptr<wGdalRaster> inras(wGdalRasterFactory::OpenFile( infile1 )) ;
	
	if( inras.get() == 0 ){
		error = infile1 + " is bad";
		return false ;
	}

	const double* inTrans  = inras->getTrans() ;//maybe not wgs84
	char* inProj = inras->getProj() ;

	double outTrans[6];
	outras->getDataset()->GetGeoTransform(outTrans) ;
	const char* outProj = outras->getDataset()->GetProjectionRef() ;


	OGRSpatialReference outSrs(outProj);
	OGRSpatialReference inSrs(inProj);

	//OGRCoordinateTransformation* pTransformer=OGRCreateCoordinateTransformation(&inSrs,
	//	&outSrs);
	OGRCoordinateTransformation* pTransformerO2I =
		OGRCreateCoordinateTransformation(
			&outSrs,
			&inSrs);

	OGRCoordinateTransformation* transFromIn2Out =
		OGRCreateCoordinateTransformation(
			&inSrs,
			&outSrs);


	const int inxsize = inras->getXSize() ;
	const int inysize = inras->getYSize() ;
	const int outxsize = outras->getWidth() ;
	const int outysize = outras->getHeight() ;

	double halfouttrans1 = outTrans[1]/2 ;
	double halfouttrans5 = outTrans[5]/2 ;
	double halfintrans1 = inTrans[1]/2 ;
	double halfintrans5 = inTrans[5]/2 ;

	int input_in_out_min_x = outxsize ;
	int input_in_out_max_x = 0 ;
	int input_in_out_min_y = outysize ;
	int input_in_out_max_y = 0 ;

	//
	cout<<"first pass..."<<endl ;
	for(int iny = 0 ; iny < inysize; ++ iny ){
		double incoory = inTrans[3] + inTrans[5] * iny + halfintrans5 ;
		for(int inx = 0 ; inx < inxsize;++inx ){
			double incoorx = inTrans[0] + inTrans[1] * inx + halfintrans1 ;

			int inval = inras->getValuei(inx,iny,0) ;
			if( inval != filldata && inval >= vmin && inval <=vmax )
			{
				double outcoorx = incoorx ;
				double outcoory = incoory ;
				transFromIn2Out->Transform(1, &outcoorx, &outcoory ) ;
				outcoorx = outcoorx - halfouttrans1;
				outcoory = outcoory - halfouttrans5;

				int outx = (outcoorx - outTrans[0])/outTrans[1] ;
				int outy = (outcoory - outTrans[3])/outTrans[5] ;

				if( outx >=0 && outx < outxsize && outy>=0 && outy<outysize )
				{
					int oldval = outras->getDatai(outx,outy,0) ;
					if( oldval == filldata ){
						outras->setDatai(outx,outy,0,inval) ;
					}else{
						//skip no min max compare, just use first value.
						// if( imethod==M_MAX) 
						// 	outras->setDatai(outx,outy,0, max( inval,oldval))  ;
						// else 
						// 	outras->setDatai(outx,outy,0, min( inval,oldval))  ;
					}
					input_in_out_min_x = min(input_in_out_min_x,outx) ;
					input_in_out_min_y = min(input_in_out_min_y,outy) ;
					input_in_out_max_x = max(input_in_out_max_x,outx) ;
					input_in_out_max_y = max(input_in_out_max_y,outy) ;
				}
			}
		}
	}
	

	cout<<"first pass done, go second pass..."<<endl ;

	int per0 =  -1;
	for(int oy = input_in_out_min_y ; oy < input_in_out_max_y ; ++ oy )
	{
		double olatcenter = outTrans[3] + outTrans[5] * oy + halfouttrans5;
		for(int ox = input_in_out_min_x ; ox < input_in_out_max_x ; ++ ox )
		{
			double oloncenter = outTrans[0] + outTrans[1] * ox + halfouttrans1 ;
			int oldval = outras->getDatai(ox,oy,0) ;
			if( oldval != filldata ) continue ;

			double incx = oloncenter ;
			double incy = olatcenter ;
			pTransformerO2I->Transform(1, &incx, &incy ) ;
			incx = incx - halfintrans1 ;
			incy = incy - halfintrans5 ;

			int iix = (incx - inTrans[0])/inTrans[1] ;
			int iiy = (incy - inTrans[3])/inTrans[5] ;

			if( iix >=0 && iix < inxsize && iiy >=0 && iiy < inysize )
			{
				double newval = inras->getValued(iix,iiy,0) ;
				if( newval != filldata && newval >= vmin && newval<=vmax )
				{
					if( oldval == filldata ){
						outras->setDatai(ox,oy,0,newval) ;
					}
					else 
					{
						//skip it , just use first value.
						// if( imethod==M_MAX) 
						// 	outras->setDatai(ox,oy,0, max( (int)newval,oldval))  ;
						// else 
						// 	outras->setDatai(ox,oy,0, min( (int)newval,oldval))  ;
					}
				}
			}			
		}
	}
	cout<<"second pass done. "<<endl ;
	OCTDestroyCoordinateTransformation(pTransformerO2I);
	OCTDestroyCoordinateTransformation(transFromIn2Out) ;
	return true ;
}

void fillGaps(wGdalRasterLiner* outras, double filldata ) {
	int filldatai = filldata ;
	const int outxsize = outras->getWidth() ;
	const int outysize = outras->getHeight() ;
	for(int oy = 1 ; oy < outysize-1 ; ++ oy )
	{
		for(int ox = 1 ; ox < outxsize-1 ; ++ ox )
		{
			int inval = outras->getDatai(ox,oy,0) ;
			if( inval == filldatai ){
				int valu = outras->getDatai(ox,oy-1,0) ;
				int vald = outras->getDatai(ox,oy+1,0) ;
				int val0 = outras->getDatai(ox-1,oy,0) ;
				int val1 = outras->getDatai(ox+1,oy,0) ;
				
				if( valu != filldatai && vald!=filldatai ){
					inval = (valu+vald)/2 ;
					outras->setDatai(ox,oy,0,inval) ;
				}else if( val0 !=filldatai && val1 != filldatai) {
					inval = (val0+val1)/2 ;
					outras->setDatai(ox,oy,0,inval) ;
				}
				
			}
		}
	}
}



int main(int argc , char** argv )
{
	cout<<"A program to composite and merge all input data into a WGS84 output data."<<endl ;
	cout<<  "v1.0.2.0 created 2022-7-14"<<endl ;
	cout<<"usage:qh5_merge infilelist.txt validmin validmax inoutfilldata out_xmin out_ymin out_xmax out_ymax out_reso out_buffer_line_count output.tif"<<endl ;
	cout<<"composite use first valid pixel, no min max compare."<<endl ;
	cout<<"infilelist.txt: one line has: one filename.tif"<<endl ;
	cout<<"e.g.:"<<endl ;
	cout<<"/some/ndvidir/ndvi1.tif"<<endl ;
	cout<<"/some/ndvidir/ndvi2.tif"<<endl ;
	cout<<"/some/ndvidir/ndvi3.tif"<<endl ;
	cout<<"the output support very large image by using wGdalRasterLiner."<<endl ;
	cout<<"only support byte,uint16,int16"<<endl ;

	GDALAllRegister() ;
	if( argc!=12 )
	{	
		cout<<"argc!=12."<<endl ;
		return 11 ;
	}

	string inlistfile = argv[1] ;
	const double vmin = atof(argv[2]) ;
	const double vmax = atof(argv[3]) ;
	const double filldata = atof(argv[4]) ;
	const double outxmin = atof(argv[5]) ;
	const double outymin = atof(argv[6]) ;
	const double outxmax = atof(argv[7]) ;
	const double outymax = atof(argv[8]) ;
	const double outreso = atof(argv[9]) ;
	int outBufferLineCount = atof(argv[10]) ;
	if( outBufferLineCount<1 ) outBufferLineCount=1 ;
	string outputfile = argv[11] ; 

	cout<<"------------"<<endl ;
	cout<<"inlistfile :"<<inlistfile<<endl ;
	cout<<"vmin :"<<vmin<<endl ;
	cout<<"vmax :"<<vmax<<endl ;
	cout<<"filldata :"<<filldata<<endl ;
	cout<<"outxmin :"<<outxmax<<endl ;
	cout<<"outymin :"<<outymin<<endl ;
	cout<<"outxmax :"<<outxmax<<endl ;
	cout<<"outymax :"<<outymax<<endl ;
	cout<<"outreso :"<<outreso<<endl ;//
	cout<<"outBufferLineCount :"<<outBufferLineCount<<endl ;//
	cout<<"outputfile :"<<outputfile<<endl ;//
	cout<<"------------"<<endl ;

	vector<string> infileVec ;
	bool listok = ReadFileList( inlistfile , infileVec ) ;
	if( listok==false ){
		cout<<"failed to read inlistfile:"<<inlistfile<<endl ;
		return 12 ;
	}

	if( infileVec.size()==0 ){
		cout<<"infileVec size equal 0. "<<endl ;
		return 13 ;
	}

	const int outxsize = (outxmax-outxmin)/outreso + 1 ;
	const int outysize = (outymax-outymin)/outreso + 1 ;
	if( outxsize<=0 || outysize<=0 ){
		cout<<"output size invalid "<<outxsize<<","<<outysize<<endl;
		return 14 ;
	}

	GDALDataType thedatatype = GDT_Byte ;
	// only one band ok.int theBandNum = 0 ;
	{
		shared_ptr<wGdalRaster> tempPtr(wGdalRasterFactory::OpenFile( infileVec[0] )) ;
		if( tempPtr.get() == 0 ){
			cout<<"bad file at 0."<<endl ;
			return 14 ;
		}
		thedatatype = tempPtr->getDataType() ;
	}

	if( thedatatype==GDT_Byte || thedatatype==GDT_Int16 
		|| thedatatype==GDT_UInt16 )
	{
		//ok
	}else{
		cout<<"not supported data type : "<<thedatatype<<endl ;
		return 15 ;
	}

	{
		wGdalRasterLiner* theOutRaster = nullptr ;
		
		wGdalRasterLinerU16 tempRasterU16 ;
		wGdalRasterLinerI16 tempRasterI16 ;
		wGdalRasterLinerByte tempRasterByte ;

		bool createok = false ;
		if( thedatatype==GDT_Byte ){
			createok = tempRasterByte.create( outputfile,outxsize,outysize,1 ) ;
			theOutRaster = &tempRasterByte;

		}else if( thedatatype==GDT_UInt16){
			createok = tempRasterU16.create( outputfile,outxsize,outysize,1 ) ;
			theOutRaster = &tempRasterU16;
		
		}else if( thedatatype==GDT_Int16){
			createok = tempRasterI16.create( outputfile,outxsize,outysize,1 ) ;
			theOutRaster = &tempRasterI16 ;

		}else{
			cout<<"not supported data type : "<<thedatatype<<endl ;
			return 16 ;
		}

		if( createok == false )
		{
			cout<<"create output raster failed."<<endl ;
			return false ;
		} 

		cout<<"filling ..."<<endl ;
		for(int iy = 0 ; iy < outysize;++iy ){
			for(int ix = 0 ; ix < outxsize;++ix ){
				theOutRaster->setDatai(ix,iy,0,filldata) ;
			}
		}
		cout<<"filling done."<<endl ;
		cout<<"set proj and transform."<<endl ;

		theOutRaster->getDataset()->SetProjection("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]"); 
		
		double outTrans[6] = {outxmin ,
			 outreso , 0.0 , 
			 outymax , 0.0 , 
			 -outreso } ;

		theOutRaster->getDataset()->SetGeoTransform(outTrans) ;

		for(int iin = 0 ; iin < infileVec.size() ; ++ iin )
		{
			cout<<"ifile "<<iin<<"/"<<infileVec.size()<<endl ;
			string error1;
			bool ok1 = FillOneFileIntoOutput(
				theOutRaster ,
				infileVec[iin] ,
				filldata ,
				vmin,vmax,
				error1) ;
			if( ok1==false ){
				cout<<"bad fill for "<<infileVec[iin]<<endl ;
				cout<<"error:"<<error1<<endl ;
				return 16 ;
			}
		}

		cout<<"all files done."<<endl ;

		//
		cout<<"last fill for gaps ..."<<endl ;
		fillGaps( theOutRaster , filldata ) ;


		cout<<"saving..."<<endl ;
	}

	cout<<"success"<<endl ;

	return 0 ;
}